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Abstract 

A newly developed method for systematically improving the convergence of path 
integrals for transition amplitudes, introduced in Phys. Rev. Lett. 94 (2005) 180403, 
Phys. Rev. B 72 (2005) 064302, Phys. Lett. A 344 (2005) 84, and expectation values, 
introduced in Phys. Lett. A 360 (2006) 217, is here applied to the efficient calculation 
of energy spectra. We show how the derived hierarchies of effective actions lead to 
substantial speedup of the standard path integral Monte Carlo evaluation of energy 
levels. The general results and the ensuing increase in efficiency of several orders 
of magnitude are shown using explicit Monte Carlo simulations of several distinct 
models. 
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1 Introduction 



Feynman's path integrals [5,6] provide the general mathematical framework for 
dealing with quantum and statistical systems. The formalism has been success- 
fully applied in generalizing the quantization procedure from the archetypical 
quantum mechanical problem of the dynamics of a single particle moving in 
one dimension, to more particles, more dimensions, as well as to more com- 
plicated objects such as fields, strings [7], etc. Symmetries of physical systems 
can be more easily treated and applied in this formalism, since it gives a simple 
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and natural setup for their use [8] . Various approximation techniques are more 
easily derived within the framework of this formalism, and it has been suc- 
cessfully used for deriving non-perturbative results. The parallel application 
of this formalism in both high energy and condensed matter physics makes it 
an important general tool [9,10]. The analytical and numerical approaches to 
path integrals have by now become central to the development of many other 
areas of physics, chemistry and materials science, as well as to the mathemat- 
ics and finance [11,12,13,14]. In particular, general numerical approaches such 
as the path integral Monte Carlo method have made possible the treatment 
of a wealth of non-trivial and previously inaccessible models. 

The key impediment to the development of the path integral formalism is a 
lack of complete understanding of the general mathematical properties of these 
objects. In numerical approaches limited analytical input generally translates 
into lower efficiency of employed algorithms. The best path generating algo- 
rithms, for example, are efficient precisely because they have built into them 
the kinematic consequences of the stochastic self-similarity of paths [15]. A 
recent series of papers [1,2,3] has for this reason focused on the dynamical im- 
plications of stochastic self-similarity by studying the relation between path 
integral discretizations of different coarseness. This has resulted in a system- 
atic analytical construction of a hierarchy of TV-fold discretized effective actions 
S$ labeled by a whole number p and built up from the naively discretized ac- 
tion in the mid-point ordering prescription (corresponding to p — 1). The level 
p effective actions lead to discretized transition amplitudes and expectation 
values differing from the continuum limit by a term of order 1/N P . 

In this paper we extend the applicability of the above method for improving 
the efficiency of path integral calculations to the evaluation of energy spectra. 
We show how the increased convergence of path integrals translates into the 
speedup in the numerical calculation of energy levels. Throughout the paper 
we present and comment on the Monte Carlo simulations conducted using the 
hierarchy of effective actions for the case of several different models including 
anharmonic oscillator, Poschl- Teller potential, and Morse potential. All the 
numerical simulations presented were done using Grid-adapted Monte Carlo 
code and were run on EGEE-II and SEE-GRID-2 infrastructure [16,17]. The 
effective actions and the codes used can be found on our web site [18]. 



2 Partition Function and Energy Spectra 

The partition function is the central object in statistical mechanics. The path 
integral formalism gives us an elegant framework for calculating partition func- 
tions which can be used either for deriving analytical approximation tech- 
niques or for carrying out numerical evaluation. The starting point is the 



2 



expression for the partition function in the coordinate basis, 



Z(0) = J daA(a,a;P) , (1) 

— oo 

where A(a, b; (5) = {b\e~@ H \a) is the quantum mechanical transition amplitude 
for going from a to b in (Euclidean) time f3. In the path integral formalism 
transition amplitudes are given as the N — > oo limit of the (TV — l)-fold integral 
expression 

1 ~ 

A N (a, b- (3) = (2^) * fdqi--- dq N ^ e" 5 ^ . (2) 



Sn is the naively discretized action of the theory, = (3/N the discrete time 
step. For the physical models that we consider the action is of the form 

p 

S = jdt (-? + V (qj) , (3) 



and its naive discretization equals 

N-l / 5 2 \ 
Sn= E [f^ + ^V(q n )) , (4) 

n=0 \ 2e JV / 



where 5 n = q n+ \ — q n , and q n = |(? n +i + Qn)- Note that we are using units 
in which the particle mass and h have been set to unity and that we are 
evaluating path integrals in the so-called mid-point ordering prescription. 

From the above we have obtained a path integral representation for the parti- 
tion function that is directly amenable to numerical evaluation. On the other 
hand, by evaluating the trace in Eq. (1) in the energy basis we find 

00 

Z((3) = e-W) = E e"^ . (5) 

?1=0 



As we can see, the partition function, or equivalently the free energy F({3), 
completely determines the energy spectrum and vice- versa. For example, if we 
define a series of auxiliary functions as 
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then it immediately follows that F n (f3) — > £?„ for large (3. It would be ideal, 
therefore, if we could calculate the free energy (and the other auxiliary func- 
tions) for arbitrarily large values of (3. This is not possible in numerical calcula- 
tions. First of all the calculations become much more demanding with growth 
of "time of propagation" (3 (just as the physics becomes more interesting). 
More importantly, when doing numerical calculations we evaluate discretized 
quantities such as Fn, and the N — > oo and (3 — > oo limits that one would 
need to perform do not commute. The best way to see this is to look at the 
free energy of an exactly solvable model - the harmonic oscillator. In this case 
the iV-fold discretized free energy (in the left ordering prescription) equals [14] 

F^(/3) = ^ In (2 sinh (£/?)), (7) 



where Cj = (2/ejv) arcsinh(u;eAr/2). This solution is illustrated in Fig. 1. It 
follows that, unlike its continuum limit F(/3), the discretized free energy F N {f3) 
does not tend to a constant value for large (3. Said another way, the discretized 
energy levels themselves depend on and thus on f3. For example, for the 
harmonic oscillator we have E Njn (ejy) = Q{n + 1/2). 
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Fig. 1. The curves depict the exact solution of the discretized free energies Fpf(f3) 
for the harmonic oscillator in the left ordering prescription given in Eq. (7) for 
various values of N. The data points give the results and error bars of the corre- 
sponding numerical calculations, used to verify the code. Parameters are u = 1 and 
N MC = 10 7 . 
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In the case of a general theory the free energy is related to its discretized 
value as F(f3) — F N (f3) = O(e^). We see that F N ((3) slowly converges to its 
continuum limit, i.e. that we need a large number of discretization points N 
to approach that value. In addition, the larger the value of (5 we want, the 
larger N must be in order to achieve a given accuracy. The price we pay is in 
the computer time which grows linearly with N . 

A recent series of papers [1,2,3] analytically studied the relation between path 
integral discretizations of different coarseness for the case of a general theory. 
This work resulted in a systematic construction of a hierarchy of iV-fold dis- 
cretized effective actions S$ labeled by a whole number p and built up from 
the naively discretized action in the mid-point prescription (corresponding to 
p = 1). The level p effective action leads to discretized transition amplitudes 
and expectation values differing from the continuum limit by a term of order 
1/N P . Thus, moving up the hierarchy we are guaranteed to get expressions 
which converge ever faster to the continuum limit. The direct application of 
these results to the free energy gives 

F(P)-F^(P) = 0(e%) . (8) 



For a given inverse temperature /3, and for < 1 the discretized free energy 
F^\f3) converges faster to the continuum as we increase the hierarchy level 
p. This is illustrated in Fig. 2. 

When using the path integral Monte Carlo method to calculate the free energy 
F(f3) there are two sources of errors. The first comes from the limited number 

— 1/2 

of Monte Carlo samples Nmc and is proportional to N M( j . The second type 
of error comes from discretization - in our case from approximating the free 
energy with F$\(3) for some N and p. As we have seen, for a given (3 this 
discretization error is proportional to N~ p . These two types of errors should 
optimally be of the same order, e.g. there is no point in decreasing the dis- 
cretization error bellow the Monte Carlo error as this would not decrease the 
overall error. In practice we fix the precision we want by choosing the number 
of Monte Carlo samples and then decrease the discretization error to match 
this either by increasing N or the hierarch level p. The second choice is far 
better; however, since computation times grow linearly with N, but are almost 
independent of p (at least for p < 9, the hierarchy levels studied in [1,2]). As 
a consequence of this, the speedup coming from using higher values of p at 
fixed precision 5 is proportional to 5~ 1+1 / p . Therefore, by using p = 9 we are 
in fact quite near to the point of optimal benefit for which the speedup of 
the new method is inversely proportional to the precision. As an illustration, 
for two decimal precision the new method gives a hundred fold speedup over 
the defining algorithm, for four decimal precision the speedup is ten thousand 
fold, etc. It is important to note that the greatest utility of the new evalua- 
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Fig. 2. The dependance of F^\(5) on iV for different levels p. The plot is for the 
anharmonic oscillator with quartic coupling g = 1, inverse temperature (3 = 1 and 
JVjuc = 10 Monte Carlo samples. The same kind of behavior is seen for other 
parameters as well as for other potentials. 

tion scheme is, therefore, when calculating quantities with high precision. We 
stress that all of this holds for < 1, i.e. as long as N > (3 is satisfied. 



3 Numerical Results 



As we have seen in the previous section, F{(3) can be evaluated with arbitrary 
precision on any interval of inverse temperatures [0, (3 max \ for any given po- 
tential by appropriately increasing and adjusting N, p, and Nmc- Let us now 
numerically compare the quality of different discretizations of the free energy 
F$ with F*, the most accurate one that may be calculated on a given set 
{Pi}. To do this we use the standard \ 2 function, 

M { ft} (AFif ) (A)) +(af*(a 



where M is the number of points in the set and AF is the Monte Carlo 
error. By including the Monte Carlo error of F* into the x 2 weights we took 
into account the fact that it is also calculated numerically. \ 2 should be around 
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one for well optimized N and p. Note that x 2 ^> 1 if the exact value of F* is 
not within the error bars of F^\ while x 2 1 if the Monte Carlo error is too 
large. 

We conducted this test on the anharmonic oscillator with quartic coupling 
V(q) = \ C L 2 + fi<7 4 - The discretized free energies were calculated for (3 G [0.5, 8] 
with step 0.5, N < 1024 and p = 1, 2, . . . , 9. The number of Monte Carlo 
samples used was 10 6 . The comparisons were done for a range of coupling 
constants g G {0,0.1,1,10,100,1000}. Taking Fy^i as the exact result, we 
calculated x 2 f° r eacn P a h °f parameters (N,p) and coupling g, and looked 
for (N,p) pairs with approximately the same values of x 2 ■ These pairs are 
given in Fig. 3. As we can see, the relation l/log 2 iV oc p that is implicit in 
Eq. (8) actually holds, i.e. the error indeed scales as N~ p . 




Fig. 3. Pairs of N and p which give similar values od x 2 ■ The plot gives l/log 2 -/V 
on y axis as a function of p. The general behavior is illustrated on the case of 
the anharmonic oscillator with quartic coupling g = 1, (3 max = 8, Nmc = 10 6 , 
X 2 ~ 2 - 4. 



We now turn to calculating the energy spectrum using the outlined efficient 
procedure for evaluating the free energy of a general theory. For the range of 
inverse temperatures j3 that will be used for numerical calculations of the 
energies we choose f3 max so that F N (j3) = F(/3) within the error bars on 
the whole [0,/3 max ] interval. We also need to ensure that all the assumptions 
mentioned above hold (ejv ^ 1, Pmax fixed). The free energy F(j3) and all its 
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auxiliary functions can be written as 
1 

1 



F n (0) = E n - i In f 1 + f) e-^~ E ^ ] . (1 ) 

i=n+l 



As a result, we have fit the numerical data to functions of the form 

F n (P) = E n -^\n(l + Ae- Bf3 ) , (11) 



where E n , A and B are the parameters of the fit. Fig. 4 shows the free en- 
ergy F(/3) (approximated by its discretization for N = 256 and p — 9) along 
with the associated auxiliary functions -F\(/3), and F 2 (/3) for the anharmonic 
oscillator with quartic coupling g — 1. Note that the class of functions given 
in Eq. (11) gives a better fit for larger values of j3. This can indeed be explic- 
itly seen from Fig. 4. The data points for the free energy F(j3) were obtained 
directly from our Monte Carlo simulations and were used to determine the 
ground state energy E . The auxiliary functions F n (j3) were obtained recur- 
sively using Eq. (6) and the already determined energy levels. The error bars 
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Fig. 4. Dependance of the free energy F and the associated auxiliary functions Fi 
and F2 on j3 for the anharmonic oscillator with quartic coupling g = 1. The solid 
lines are the fits to curves of the form given in Eq. (11). The horizontal lines in 
black correspond to the energy levels E n determined from these fits (see Table 1). 
Numerical simulations were performed with p = 9 level improved actions, N = 256, 
and N M c = 10 7 . 
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presented in the figure also follow directly from Eq. (6) and are given by 



For large inverse temperatures (3 the above denominator becomes exponen- 
tially small, and so the error bars become very large. Such points soon cease 
to give relevant contributions to the calculations of the corresponding energy 
level owing to the fact that we use a weighted fit. Note that, in fact, the lack 
of exponential growth of error bars with (5 is an indication of bad data points! 

This effect of growing error bars becomes more pronounced for higher energy 
levels. In addition, from Eq. (12) we see that there is an accumulation of 
errors associated with all the lower energy levels. Both of these effects taken 
together give practical limits to the number of energy levels we can calculate. 
The precise depth to which we can probe the energy spectrum depends on the 
number of Monte Carlo samples used as well as the number of points $ selected 
within the range of inverse temperatures available to us. As an illustration, 
Table 1 gives the low lying energy levels of the anharmonic oscillator for several 
values of coupling g. For all of these calculations we use the same range of (3. 
The ground state energy level was calculated to five significant digits for all 
values of g. As we have already noted the errors increase as we go to higher 
energy levels. In fact, this increase is faster for larger couplings since then 
the energies themselves become higher and so the e^^ En terms become much 
smaller. 
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1000 


2.3578(2) 









Table 1 

Low lying energy levels of the anharmonic oscillator with quartic coupling g, calcu- 
lated using N = 256, p = 9, and N MC = 10 7 . 



We have conducted explicit Monte Carlo calculations of the spectra of the 
Poschl- Teller and Morse potentials and have obtained the same qualitative be- 
havior. In particular, we have explicitly determined that the expected speedup 
in convergence, coming from using the p-level hierarchy of effective actions, 
holds for all of these potentials. 
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Table 2 

Low lying energy levels of the modified Poschl- Teller potential, calculated using 
N = 256, p = 9, and N M c = W 7 . 

Obtained low lying energy levels for several values of the parameters of the 
modified Poschl- Teller potential, 

^ 2 cosh 2 ax ' ^ 



are given in Table 2. We considered this exactly solvable potential since it 
allows comparison of numerically calculated energy levels and the exact ones, 
given by 

a 2 

E exact = (\ - l-n) 2 ,0<n<A-l,neN. 

As can be seen from Table 2, numerical results are in excellent agreement with 
the exact energy levels even for a small value of discretization coarseness N. 

As a conclusion, we have investigated a newly developed method for increasing 
the convergence of path integrals to the continuum limit. The method has 
previously been shown to lead to a many order of magnitude speedup in the 
numerical evaluation of path integrals for transition amplitudes [1,2,3] and 
expectation values [4]. In this paper we have applied that method to the 
evaluation of energy spectra. We have shown that the above stated increase 
in convergence leads to a significant increase of the efficiency of path integral 
Monte Carlo calculations of low lying energy levels of a generic theory. The 
analytical results were checked explicitly in a series of Monte Carlo simulations 
of several distinct models over a wide range of parameters. 
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